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Abstract 

A highly efficient formulation of moment equations for stochastic reaction networks is introduced. 
It is based on a set of binomial moments that capture the combinatorics of the reaction processes. 
The resulting set of equations can be easily truncated to include moments up to any desired 
order. The number of equations is dramatically reduced compared to the master equation. This 
formulation enables the simulation of complex reaction networks, involving a large number of 
reactive species much beyond the feasibility limit of any existing method. It provides an equation- 
based paradigm to the analysis of stochastic networks, complementing the commonly used Monte 
Carlo simulations. 

PACS numbers: 05.10.Gg, 82.40.-g, 82.40. Qt, 02.50.Fz, 02.50. Ga 



Stochastic reaction networks appear in many natural systems, from grain-surface chem- 
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3| and to ecological systems 
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istry in the interstellar medium to biochemistry in cells 
^ . In order to characterize these networks one wishes to evaluate the time dependent pop- 
ulations of the reactive species. Due to the large fluctuations in stochastic systems, the 
rate equations fail and stochastic methods 5|, |6[ are required. These methods are based on 
the master equation, which can be solved either by direct integration or by Monte Carlo 
simulations {7,8]. The problem with these methods is that they quickly become infeasible as 
the number of reactive species increases. More specifically, the number of equations in the 
master equation proliferates exponentially with the number of species, making the stochastic 
analysis difficult even for empirical networks of moderate size 9 
is to derive moment equations by tracing over the master equation 
one directly computes different moments related to the populations of the reactive species. 
The problem is that higher order moments appear on the right hand side of the equations, 
making it difficult to obtain a closed set of equations. In most closure schemes, only the 
first and second moments are included 12|. 

In this Letter we present a general and highly effective formulation of moment equations, 
expressed in terms of the binomial moments, defined below. These are linear combinations 
of the ordinary moments, that capture the structure of the reactions and provide much 
physical insight. As a result, this formulation enables the inclusion of moments up to any 
desired order. The equations are linear and their number is dramatically reduced compared 
to the master equation, enabling the analysis of complex reaction networks. 

Consider a reaction network, consisting of J reactive species, Xi, i = 1, . . . , J. These 
species undergo reactions of the form 
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J^n^Xi ^ J]m^Xj, (1) 

i=i j=i 

where the stoichiometric coefficients rij and mj are integers. The order of the reaction is 
given hj n = Yli=i^iy namely the number of reactants on the left hand side of Eq. ([1]). In 
the context of chemical reactions it is common to omit reactions of order higher than n = 2. 
However, in the formulation presented below there is no such limitation. The reactions 
presented in Eq. ([T]) can be expressed in a vector form as n — )■ m, namely a set of molecules 
with stoichiometric coefficients n, reacts to form a new set m. The rate constant for the 
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molecular configuration n to react is given by Tn (s~^). Certain reactions may take several 
paths, with different probabilities or branching ratios. The probability that the reaction will 
result in the configuration m is denoted by P^. These probabilities satisfy J2m^H ~ ^■ 
Thus, the rate constant for the reaction n — )■ m is TfiP^^. 

In the stochastic analysis one describes the state of the system by the vector N = 
{Ni,...,Nj), where Ni is the number of copies of species Xi. Below we introduce a 
combinatorial approach, which will later allow us to express the master equation and 
the moment equations in a transparent form. Let v = {vi, . . . ,vj) be a vector of inte- 
gers. It can be expressed as a linear combination of the basis vectors, ei,...,ej, where 
Cj = (0, . . . , = 1, . . . , 0). We denote by Q^j a combination of molecules consisting of ex- 
actly Vi copies of the species Xi. The number of such combinations that exist in a system 
at the state N is given by 

W{N,v)=(^^_)j, (2) 

where (^) = nf=i (C^') CH') binomial coefficient. To illustrate the motivation for 

this definition, consider the reaction n ^ m. It occurs at a rate proportional to Tf^P^, and 
to the number of Qa combinations which are present in the system, given by W(A^, n). 

— # — * 

Let P{N) represent the time dependent probability for the system to be in the state N. 
The master equation for P{N) takes the form 



dP{N) 



W{N + n-m,n)P{N + n-m) -^JV{N,n)P{N) . (3) 
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In this equation one sums over all the reactions n — )■ m in the system. These reactions 
yield a positive contribution to P{N) when the state of the system is -|- n — m, and a 
negative contribution to P{N) when the system is in the state N. In numerical simulations 
the master equation must be truncated in order to maintain a finite number of equations. 
This is achieved by setting upper cutoffs Cj, z = 1, . . . , J, such that P{N) = if Ni > Ci for 
any value of i. The number of coupled equations, A''^; = n/=i(^* + 1)' grows exponentially 
with the number of reactive species, J. This severely limits the applicability of the master 
equation to complex reaction networks jol, lo| . 

A more compact description of stochastic reaction networks can be obtained using mo- 
ment equations. These equations, derived by tracing over the master equation, consist of 



ordinary differential equations for the time derivatives of the moments (A'^"^ ■ ■ ■ N^), where 
tti are integers. The order of a specific moment is given hj k = Yli=i The difficulty with 
moment equations is that higher order moments appear on the right hand side of each equa- 
tion. To obtain a closed set of equations one needs to apply a suitable truncation scheme, 



in which the higher order moments are expressed in terms of low order moments |ll|-|l4l|. In 
practice, the truncation is typically done at the level of third order moments, namely only 
the first and second order moments are taken into account. Making the truncation at higher 
orders turns out to be prohibitively complicated even for relatively simple networks. Below 
we introduce a different formulation of the moment equations, which enables to extend the 
truncation up to any desired order. 

Consider a reaction network described by the probability distribution P{N). The bino- 
mial moment (W^) is defined as the average number of combinations of the form that 
appear in the system. It is given by 

{W^) = J2w{N,v)P{N). (4) 

TV 

To understand the meaning of the binomial moments, consider the case where u = Cj. Here 
the corresponding binomial moment, (W^.), is given by the average number of combinations 
of a single copy of the species Xj. This is simply the average population size, (Ni). In case 
that V = Ci + ej, the corresponding moment is (W^;) = (NiNj), which stands for the average 
number of Xi-Xj pairs present in the system. Note that for v = 2ej the corresponding 
binomial moment is (W^) = {{Nf) — {Ni))/2. Similarly, the number of Xi triplets is given 
by (W^) = ((iVf) -3(iVf) + 2(iVi))/6, where v = Se^. The order the binomial moment (W^) 

To obtain the binomial moment equations we express the time derivative of (W^) using 
Eq. (jll), where the time derivative of P{N) is taken from the master equation [Eq. ([3])], and 



a summation is taken over 15| . The resulting binomial moment equations take the form 



dt ^ 



(8) 



where 
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The first term in Eq. ([5]) accounts for positive contributions of the reactions to (W^), 
while the second term acounts for the negative contributions. We first consider the positive 
contributions. Consider a single combination Qt7+^„^. There are on average (W^+^_^) 
such combinations in the system. For each one of these combinations, the reaction n ^ rh, 
produces a new combination, so (W^g) increases. The rate constant accounts for the 
rate of formation of combinations by the reactions n ^ w for all possible choices of w. 
The reaction rate is also proportional to the binomial coefficient (^^'^~'") , that accounts for 
the number of combinations Qh in each of which may undergo the reaction. 

We now refer to the negative contributions. Consider the combination Q^+rn which 
undergoes a reaction of the form n ^ for any possible choice of w. The overall rate of these 
reactions is given by Tfii^'^S^^ . Each time such a reaction takes place, a single combination 
is removed from the system. The removed combination can be decomposed into and 
Qn-rh- Note that there are (^) different possibilities to perform this decomposition of Q^. 
When QtJi is removed, the combination Q^+^n is replaced by by Q^. This combination is 
then eliminated by the subsequent removal of Qn-m- 

Eq. ([5]) is not in a closed form, because higher order moments appear on the right hand 
side of each equation. However, the binomial moments tend to decrease as their order 
increases. To demonstrate this feature, consider a binomial moment (W^;) of order k = 
Yli=i '^i- It represents the average number of appearances of a certain combination Q^j which 
consists of k molecules in the system. In a small system, where the average copy numbers are 
small, (Wi/) tends to decrease as k increases. This enables us to use the following truncation 
scheme. We choose a cutoff C such that (W^) is set to zero whenever k > C. The number 
of different binomial moments of order k is given by C^^^^^)- Thus the number of binomial 
moment equations, after the truncation is carried out, is Ne = ^^Li {''~^k~^)' While the 
number of equations in the master equation grows exponentially with J, the number of 
binomial moment equations scales only polynomially with J. Moreover, in practice one can 
obtain accurate results using surprisingly low values of the cutoffs. In fact, in an earlier 
version of this formulation, it was shown that for a broad range of conditions the cutoff 
C = 2 is sufficient [isl, In this case, the equations include first order moments that 
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account for the average population sizes and second order moments that account for the 
number of pairs of species Xi and Xj in the system, from which the reaction rates are 
evaluated. In cases where a cutoff of C = 2 is not sufficient, one may raise the cutoff until 
accurate results are obtained. In the small system limit the required cutoff is usually low. 
In the large system limit, stochastic equations are no longer required and can be replaced 
by the rate equations. 

To demonstrate the method, we apply it to the reaction network shown in Fig. [1] This 
network consists of processes involving single molecules as well as pairs of molecules. The 
network includes 10 reactive species, 3 zero order reactions, 14 first order reactions and 12 
second order reactions. The zero order reactions lead to the formation of Xi, X^ and X3 
molecules, where Pq' = 1/3 for i = 1,2 and 3. The rest of the species are formed via first 
and second order reactions. The first order reactions include the degradation of each of the 
reactive species, and the dissociation of Xq and X7. The four first order reactions involving 
Xq are given by ee — )■ (degradation), cq — )■ 2ei + 64, ee — )■ 64 + eV and ^ ei + e^, where 
P5 = 0.999 and p|"^i+^^* = p|'*+'^7 = pf^^'^ jl = 0.00025. The two first order reactions 

eg 66 £6 £6 ' 

involving X^ are eV — and eV 2ei, where P5, = 0.999, P^^^ = 0.001. The second order 
reactions involve all the pairs of nodes connected by edges. The reaction of Xi and X4 
includes three reaction paths, ei + 64 — )■ 65, ei + 64 — 63 + eV and 61 + 64—7' 2ei + 63, where 
P5\. = 0.25, P5^+5 = 0.5 and P^'it'^^ ^ o.25. The paths for the reaction of X5 and Xq 

£1+64 ' 61+64 ei+64 ^ o u 

are 65 + ce — )■ 5ei + 2e3, 65 + ce — )■ 65 + and 65 + ce — )■ eg, with the probabilities 1/4, 1/4 
and 1/2, respectively. This means that upon encounter of a pair of X5 and Xq molecules, 
they either dissociate into their fundamental components, remain unchanged or combine to 
form the molecule Xg. To characterize the size of the system, we introduce the parameter 
5*. The rate of zero order reactions is taken to be proportional to the size of the system, 
while the rate of second order reactions is inversely proportional to the size of the system 
[l^. We thus set Tq = gS, T^. = di and Tg^+ej = clij/S. The parameters we use are g = 1 
s~^ for the zero order reactions, di = 0.3 s~^, i = 1, . . . , 3, for the first order reactions and 
Qij = 1 for the second order reactions that are included in the network. 

For this reaction network, setting a cutoff of C = 3 for all the species, one obtains 
(3 + 1)^*^ ~ 10^ equations in the master equation. In contrast, the same cutoff set in the 
binomial moment equations yields only 285 equations. A lower cutoff of C = 2 results in 65 
equations, compared with approximately 6 x 10^ equations in the master equation. After 
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solving the binomial moment equations, one can extract the average population sizes, which 
are given by the (W^.) moments. 

In Fig. [2] we show the population sizes of several selected species vs. the system size, S. 
These results were obtained from the binomial moment equations with C = 2 (circles), under 
steady state conditions. For small systems, the results are in perfect agreement with those 
of the master equation (solid lines). The master equation results were obtained using the 
Gillespie algorithm 0], since direct integration of the master equation is already infeasible 



in this case 



171]. For large systems, the stochastic results converge to those obtained from 



the rate equations (dashed lines). For most species, the results of the binomial moment 
equations with such low cutoff, coincide with those of the master equation not only in the 
small system limit but also for S ^ 1. In fact, for certain species, it is sufficient to choose 
a cutoff of C = 2 to account for the abundances in the entire range of system sizes. This 
surprisin g fea ture was discussed in an earlier formulation of the moment equations, presented 
in Refs. [isl . [l^ . In case that the the cutoff of C = 2 in insufficient one may raise it to 
3 or 4, until a smooth convergence to the rate equation results is obtained in the large 
system limit. Results obtained for a cutoff of C = 3 (+) are shown for (Ni) and (A^2)- 
Consider the population of the species Xg in the system. This species is produced by the 
reaction 65 + 66 — )■ eg. Thus, in order for Xg to be produced, a pair of and Xq must 
be simultaneously present in the system. However, such pairs form only when there are 
triplets in the system. For instance, when the combination Q^^+^^+^g transforms into Qg^+^Q 
through the second order reaction 61 + 64 — 65. Thus, a cutoff lower than C = 3 in the 
binomial moment equations would terminate the production of Xg even in the stochastic 
limit. In such cases, the ability to extend the equations to higher cutoffs is crucial. In Fig. 
[3] we show the convergence of the binomial moment equations to the rate equation results 
(dashed lines) for the average population of X2. Results are shown for cutoffs of C = 2 
(circles), C = 3 (+), C = 4 (triangles), C = 5 (x) and C = 6 (squares). Already for C = 5 
or 6 the convergence to the deterministic results is smooth. When a higher cutoff is required, 
one can safely use the rate equations. 

In conclusion, the binomial moment equations provide a highly efficient equation-based 
methodology for the simulation of stochastic reaction networks. The binomial moments 
capture the essence of the combinatorics that governs the reaction rates, reflecting the stoi- 
chiometric structure of the reactions. The closure scheme is fully controlled and determined 
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by the maximal number of particles allowed to reside simultaneously in the system, which 
is clearly limited by the system size. Unlike the ordinary moment equations, the binomial 
moment equations can be constructed to any order using an automated procedure, taking 
the network structure as input. The number of equations is dramatically reduced. This 
method opens the way to systematic studies of large and complex stochastic networks be- 
yond the feasibility limit of existing methods. Moreover, as an equation-based paradigm, it 
is amenable to analytical treatments that are expected to provide crucial insight about the 
networks. 
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FIG. 1: (Color online) A reaction network consisting of 10 reactive species. First order reactions 
appear as arrows connecting the reacting species to its products. Second order reactions are 
denoted by edges connecting the two reacting species. Arrows are drawn from the edge to the 
reaction products. Different reaction paths are marked by the indices appearing by the arrows. 
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FIG. 2: (Color online) The average populations sizes, (Ni), of several species from the reaction 
network appearing in Fig. [H vs. the system size, S, as obtained from the moment equations under 
steady state conditions (circles). These results were obtained for a cutoff of C = 2. In the limit of 
small system sizes the results are in perfect agreement with the results obtained from the master 
equation (solid lines). The rate equation results (dashed lines) show significant deviations for small 
sizes. For the species Xi and X2 the moment equations deviate for large systems, and in order to 
get accurate results one has to use a higher cutoff. For these species we display results obtained 
from the moment equations with a cutoff of C = 3 (+). The species Xg cannot be produced without 
the presence of triplets in the system (bottom right display). A cutoff of 2 is thus insufficient for 
obtaining its population size. However using the moment equations with a cutoff of 3 (+) provides 
accurate results for the entire range of system sizes. 
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FIG. 3: (Color online) Here we focus on the transition between the stochastic regime S < 1 and 
the deterministic regime S > 1. Results are shown for obtained from the moment equations 
with a cutoff of 2 (circles), 3 (+), 4 (triangles), 5 (x) and 6 (squares). As the cutoff used in the 
moment equations is raised, the convergence to the rate equation results (dashed line) becomes 
smoother. When a cutoff higher than 6 is required, one reliably enters the range of validity of the 
rate equations. 
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